A NONLINEAR ELASTICITY MODEL OF MACROMOLECULAR 
CONFORMATIONAL CHANGE INDUCED BY ELECTROSTATIC FORCES 
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Abstract. In this paper we propose a nonlinear elasticity model of macromolecular 
conformational change (deformation) induced by electrostatic forces generated by an 
implicit solvation model. The Poisson-Boltzmann equation for the electrostatic poten- 
tial is analyzed in a domain varying with the elastic deformation of molecules, and a new 
continuous model of the electrostatic forces is developed to ensure solvability of the non- 
linear elasticity equations. We derive the estimates of electrostatic forces corresponding 
to four types of perturbations to an electrostatic potential field, and establish the exis- 
tance of an equilibrium configuration using a fixed-point argument, under the assumption 
that the change in the ionic strength and charges due to the additional molecules causing 
the deformation are sufficiently small. The results are valid for elastic models with arbi- 
trarily complex dielectric interfaces and cavities, and can be generalized to large elastic 
deformation caused by high ionic strength, large charges, and strong external fields by 
using continuation methods. 
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1. An Electro-Elastic Model of Conformational Change 

A number of fundamental biological processes rely on the conformational change of 
biomolecules and their assemblies. For instance, proteins may change their configu- 
rations in order to undertake new functions, and molecules may not bind or optimally 
bind to each other to form new functional assemblies without appropriate conformational 
change at their interfaces or other spots away from binding sites. An understanding of 
mechanisms involved in biomolecular conformational changes is therefore essential to 
study structures, functions and their relations of macromolecules. Molecular dynamics 
(MD) simulations have proven to be very useful in reproducing the dynamics of atomistic 
scale by tracing the trajectory of each atom in the system [37]. Despite the rapid progress 
made in the past decade mainly due the explosion of computer power and parallel com- 
puting, it remains a significant challenge for MD to study large-scale conformational 
changes occurring on time-scales beyond a microsecond [5]. Various coarse-grained 
models and continuum mechanics models are developed in this perspective to comple- 
ment the MD simulations and to provide computational tools that are not only able to 
capture characteristics of the specific system, but also able to treat large length and time 
scales. The prime coarse-grained approaches are the elastic network models, which de- 
scribe the biomolecules to be beads, rods or domains connected by springs or hinges 
according to the pre-analysis of their rigidity and the connectivity. Elastic network mod- 
els are usually combined with normal mode analysis (NMA) to extract the dominant 
modes of motions, and these modes are then used to explore the structural dynamics at 
reduced cost [9]. Continuum models do not depend on these rigidity or connectivity 
analysis. On the contrary, the rigidity of the structure shall be able to be derived from the 
results of the continuum simulations. Typical continuum models for biomolecular sim- 
ulations include the elastic deformation of lipid bilayer membranes [35] and the gating 
of mechanosensitive ion channels [34] induced by given external mechanical loads. It is 
expected that with more comprehensive continuum models we will be able to simulate 
the variation of the mechanical loads on the macromolecules with their conformational 
change, and investigate the dynamics of molecules by coupling the loads and deforma- 
tion. This article takes an important step in this direction by describing and analyzing the 
first mathematical model for the interaction between the nonlinear elastic deformation 
and the electrostatic potential field of macromolecules. Such coupled nonlinear models 
have tremendous potential in the study of configuration changes and structural stability 
of large macromolecules such as nucleic acids, ribosomes or microtubules during various 
electrostatic interactions. 

Our model is described below. Let G be a smooth open domain whose boundary 
is noted as dVl; see Fig.(l). Let the space occupied by the flexible molecules ^Imj be a 
smooth subdomain of fi, while the space occupied by the rigid molecule(s) is denoted by 
flmr- Let the remaining space occupied by the aqueous solvent be ^Is- The boundaries 
of ilmf and ^Imr are denoted by Fj and F^, respectively. We assume that the distance 
between molecular surfaces and 

mm{\x - y\ : X eTfUTr,y e dfl} (1.1) 

is sufficiently large so that the Debye-Huckel approximation can be employed to de- 
termine a highly accurate approximate boundary condition for the Poisson-Boltzmann 
equation. There are charges atoms located inside Vlmf and ^Imr, and changed mobile ions 
in fis- The electrostatic potential field generated by these charges induces electrostatic 
forces on the molecules ^Imf and ^Imr- These forces will in turn cause the configuration 
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change of the molecules. We shall model this configuration rearrangement as an elas- 
tic deformation in this study. Specifically, we will investigate the elastic deformation 
of molecule ^Imf (which is originally in a free state and not subject to any net external 
force) induced by adding molecule and changing mobile charge density in rig. This 
body deformation leads to the displacement of charges in Vlmf and the dielectric bound- 
aries, which simultaneously lead to change of the entire electrostatic potential field. It is 
therefore interesting to investigate if the deformable molecule ^Imf has a final stable con- 
figuration in response to the appearance of fi^r and the change of mobile charge density. 




Figure 1 . Illustration of macromolecules immersed in aqueous solvent environment. 

Within the framework of an implicit solvent model which treats the aqueous solvent in 
fls as a structure-less dielectric, the electrostatic potential field of the system is described 
by the Poisson-Boltzmann equation (PBE) 

Nf+Nr 

- V ■ (eV0) + K^sinh(0) = qiS{xi) in Q, (1.2) 

i 

where 5{xi) is the Dirac distribution function at Xi, Nj + Nr is the number of singular 
charges of the system including the charges in Vlmf (i.e. Nj) and Vlmr (i-e-, Nr). The di- 
electric constant e and the modified Debye-Hiickel parameter k are piecewise constants 
in domains fim/, ^mr and i^^. In particular, k = in and i^^r because it models the 
free mobile ions which appear only in the solvent region Qg. The dielectric constant in 
the molecule and that in the solvent are denoted with and e^, respectively. Readers are 
referred to [26, 27] for the importance of the Poisson-Boltzmann equation in biomolecu- 
lar electrostatic interactions, and to [2, 3, 4, 28, 29, 30, 31] for the mathematical analysis 
as well as various numerical methods for the Poisson-Boltzmann equation. 

The finite(large) deformation of molecules is essential to our coupled model, but can 
not be described by a linear elasticity theory. We therefore describe the displacement 
vector field u(x) of the flexible molecule Vlmf with a nonlinear elasticity model: 

-div(T(u)) = f, in 

T(u) ■ n = f, on rj 

where f;, is the body force, fs is the surface force and T(u) is the second Piola-Kirchhoff 
stress tensor. In this study we assume the macromolecule is a continuum medium obey- 
ing the St Venant-Kirchhoff law, and hence its stress tensor is given by the linear(Hooke's 
law) stress-strain relation for an isotropic homogeneous medium: 

T(u) = (I + Vu)[ATr(E(u))I + 2^E(u)]. 
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Here A > and > are the Lame constants of the medium, and 

E(u) = ^(Vu^ + Vu + Vu^Vu) 

is the nonlinear strain tensor. The equation (1.3) is nonlinear due to the Piola transfor- 
mation (I + Vu) in T(u), and the quadratic term in the nonlinear strain E(u). The third 
potential source of nonlinearity, namely a nonlinear stress-strain relation, is not consid- 
ered here; however, our methods apply to this case as well. 

It is noted that Eq.(1.3) is defined in the undeformed molecule body il'^j with unde- 
formed boundary r°, while the Poisson-Boltzmann equation (1.2) holds true for real 
deformed configurations. The deformed configuration is unknown before we solved 
the coupled system. We therefore define a displacement-dependent mapping ^(u){x) : 
— 7- and apply this mapping to the the Poisson-Boltzmann equation such that it can 
also be analyzed on the undeformed molecular configuration. In this map $(u)(x) 
is X + u where X is the identity mapping. A key technical tool in our work is that this 
mapping is then harmonically extended to ^7 to obtain the maximum smoothness. Apply 
this mapping, the Poisson-Boltzmann equation 1 .2 changes to be 

Nf+Nr 

- V ■ (eF(u)V0) + J(u)/s:2 sinh(0) = ^ J{\x)q^{^{x) - ^{x,)) in fi, (1.4) 

i 

where J(u) is the Jacobian of $(u) and 

F(u) = (V<l>(u))-V(u)(V$(u))-^. (1.5) 

This matrix F is well defined whenever $(u) is a C^-diffeomorphism [6]. The functions 
in Eq.(1.4) should be interpreted as the compositions of respective functions in Eq.(1.2) 
with mapping $(a;), i.e., 0(x) = 0($(a;)), e{x) = e($(x)) and k{x) = k($(x)). 

In this paper, we shall analyze the existence of the coupled solution of the elastic- 
ity equation (1.3) and the transformed Poisson-Boltzmann equation (1.4). These two 
equations are coupled through displacement mapping $(u) in the Poisson-Boltzmann 
equation and the electrostatic forces to be defined later. The solution of this coupled 
system represents the equilibrium between the elastic stress of the biomolecule and the 
electrostatic forces to which the biomolecule is subjected. The existence, the unique- 
ness and the 11/^'^ -regularity of the elasticity solution have already been established by 
Grandmont [6] in studying the coupling of elastic deformation and the Navier-Stokes 
equations; thus in this work we shall focus on the solution to the transformed Poisson- 
Boltzmann equation and to the coupled system. We shall define a mapping S from an 
appropriate space Xp of displacement field u into itself, and seek the fixed-point of this 
map. This fixed-point, if it exists, will be the solution of the coupled system. A critical 
step in defining 5* is the harmonic extension of the Piola transformation from ^Imj to 
f2 and M'^. The regularity of the Piola transformation determines not only the existence 
of the solution to the transformed Poisson-Boltzmann equation, but also the existence 
of the solution to the coupled system. Because most of our analysis will be carried out 
on the undeformed configuration we will still use i^mf, ^mr, ^s, T/, to denote the un- 
deformed configurations of molecules, the solvent and the molecular interfaces, unless 
otherwise specified. 

The paper is organized as follows. In Section 2 we review a fundamental result 
concerning the piecewise VT^'^-regularity of the solutions to elliptic equations in non- 
divergence form and with discontinuous coefficients. The nonlinear elasticity equation 
will be discussed in Section 3, where the major results from [6] are presented without 
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proof. The Piola transformation will be defined, harmonically extended, and then an- 
alyzed. In Section 4 we will prove the existence and uniqueness of the solution to the 
Piola-transformed Poisson-Boltzmann equation, generalizing the results in [2] for the 
un-transformed case. Both L°° and W'^'''^ estimates will be given for the electrostatic 
potential in the solvent region, again generalizing results in [2]. We will then define the 
electrostatic forces and estimate these forces by decomposing them into components cor- 
responding to four independent perturbation steps. The estimates of these components 
are obtained separately and the final estimate of the surface force is assembled from these 
individual estimates. The coupled system will be finally considered in Section 6 where 
the mapping S will be defined, and the main result of the paper will be established by 
applying a fixed-point theorem on this map to give the existence of a solution of the 
coupled system. 

2. Notation and Some Basic Estimates 

In what follows W'''^{'D) will denote the standard Sobolev space on an open domain 
V, where V can be f2, ilm or rig. While solutions of the Poisson-Boltzmann have low 
global regularity in we will need to explore and exploit the optimal regularity of the 
solution in any subdomain of Vl. For this purpose, we define W^'^(^7) = W'^'^{^lm) + 
W'^'''{fls) where + is the direct sum. Every function G W^'^ can be written as = 
4>m{x) + 4>s{x) where (f)m{x) G W^^'^(fim), 4>s{x) G W'^'''{Qs), and has a norm 

||0||w2,p = ||0m||w2,p(Q^) + ||0s||v^2,p(Q^). (2.1) 

Similarly, we define a class of functions C = C{^1) which are continuous in either sub- 
domain and may have finite jump on the interface, i.e., a function a G C is given by 
a = am + CLs where G C(f2m), «s ^ C{ils) are continuous functions in their respec- 
tive domains. The norm in C is defined by 

llc^llc = ||am||c(nm) + ||'^s||c(f7s)- 

We recall two important results. The first is a technical lemma which will be used for 
the estimation of the product of two W^''' functions; this is sometimes called the Banach 
algebra property. 

Lemma 2.1. Let 3<p<oo,l<q<pbe two real numbers. Let Q be a domain in M.^. 
Let u G W^''p{VL),v G W^^'^iyL), then their product uv belongs to W^''^, and there exists 
a constant C such that 

\\uv\\w^.i{n) < C'||M||vi^i,P(n)||'i^||vi/i.9(n). 

For the proof of this lemma we refer to [1]. In this paper we will apply Lemma (2.1) 
to the case with p = q. The second result is a theorem concerning the estimate of 
elliptic equations with discontinuous coefficients. 

Theorem 2.2. Let Q and Qi GG Q be bounded domains ofM.^ with smooth boundaries 
dQandT. LetQi = (QiUT) andVL2 = Let A be a second order elliptic operator 

such that 
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Then there exists a unique solution u G W^'^ for the interface problem 

{Au){x) = finVL 
[u] = U2 — Ui = on r 
[Bun] = B2VU2 ■ n — BiVui ■ n = h onV 

u = g on dQ 

providing that aik G C{n),Bi G C{T)J e LP{n),ge W^-^/P'P{dn),h e W^-^/P'P{T), 
where n is the outside normal to Qi. Moreover, the following estimate holds true 

\\u\\w'2.p{n) < K (||/||LP(n) + ||^||vyi-i/p.p(r) + \\g\\w^-^/P'P{dn) + ||'"IUp(!^)) , (2.2) 
where the constant K depends only on fi, ^2, p and the modulus of continuity of A. 

Theorem (2.2) is fundamental to various results about elliptic equations with discon- 
tinuous coefficients; For example, the global regularity and estimates of Babiiska 
[20], the finite element approximation of Chen et al. [21], a prior estimates for second- 
order elliptic interface problems [22], the solution theory and estimates for the nonlinear 
Poisson-Boltzmann equation [2], and the continuous and discrete a priori L°° estimates 
for the Poisson-Boltzmann equation along with a quasi-optimal a priori error estimate 
for Galerkin methods [2] applied to the Poisson-Boltzmann equation. For the proof of 
Theorem (2.2) and the more general conclusions for high-order elliptic equations with 
high-order interface conditions we refer to [23, 24]. 

3. Nonlinear Elasticity and the Piola Transformation 

We first state a theorem concerning the existence, uniqueness, regularity and the esti- 
mation of the solution to the nonlinear elasticity equation [6] : 

Theorem 3.1. Let the body force tb G L'^iVLmf) and the surface force is ^ H^^"^/^'^(r/), 
where 3 < p < 00. There exists a neighborhood of in L^^Qraf) x W^'^^^'^iT f) 
such that if{fb, fs) belongs to this neighborhood then there exists a unique solution u G 

w''P{nmf)nw^;^jn^j)of 

-div(T(u)) = ffe in Vlmf, 

T{u)n = isOn r/\r/o, 

u = on Tfo, (3.1) 

/ (1 + Vu)J(u)(I + Vu)-^-n = 3|l]™/|, 

where T jq is a subset ofTf equipped with homogeneous Dirichlet boundary condition, 
1 is the unit matrix. The last equation represents the incompressibility condition of the 
elastic deformation. Moreover, the solution can be estimated with respect to the force 
data: 

\\u\\w^,p(Q.^i) < C{\\fb\\LP{n^f) + \%\\w^-i/P,P(rf))- (3-2) 
Proof See [7]. □ □ 

Remark 3.2. It is noticed that u G C^'^~^/P(f2m/) because of the continuous embedding 

ofW^PiQ^f) in C'''-yp{Umf)forp > 3. 

The displacement field u{x) solved from Eq.(3.1) naturally defines a mapping $(u) = 
X-l-u in Qjnf where X is the identity mapping. This mapping $(u)(x) has to be appropri- 
ately extended into M'^ \ ^l^f to yield a global transformation for the Poisson-Boltzmann 
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equation. It is critical in what follows that this extension has various favorable properties, 
which leads us to define a global mapping by harmonic extension: 

$(u) = /^ + " "^^-^ (3.3) 
1 X + w otherwise 

where w solves 

Aw = Oin W^XUmf, 

^ ^ (3.4) 
w = u on Fj. 

The following crucial lemma concerns the regularity of $(u) and the invertibility of 
V$(u): 

Lemma 3.3. Let $(u) be defined in Eq.(3.3), we have 

(a) <l>(u) G W^'Pi^mf) and <l>(u) G C^{R^ \ H^/). 

(b) There exists a constant M > such that for all ||u||^2,p(q) < M, V$(u) is an 

invertible matrix in W^'P{Qrnf) and in C°^(]R^ \ ^Imf)- 

(c) Under condition (b), ^{x) is one-to-one on M^, and is a C^-diffeomorphismfrom 
^Imf to <l>(u)(f2m/), and is a C°°-diffeomorphismfrom \ ilmf to $(u)(M^ \ 

Proof. That $(u) G W'^'P{nmf) follows directly from its definition. Also, of <l>(u) G 
C°°(M^ \ ^mf) since $(u) = X + w while w is harmonic hence analytical in $(u) G 
(joo (]^3^^^^^ because it is the solution of the Laplace equation (3.4). For the invertibility 
of $(u) in W^'P{rirnf) we refer to Lemma 2 in [6] or Theorem 5.5.1 in [7], which says 
that if a u G fi^/ is differentiable and 

|Vu(x)| < C 

for some constant depending on ^Imf, then V$(u) = I + Vu > V x G ^mf and 
I + Vu is injective on ^Imf- The invertibility of V$(u) therefore follows from the facts 
that u G C'^''^'^/P(nmf) such that for sufficiently small M 

|Vu| < \\u\\ci,i-3/p(^n^f) < Ci||u||vt/2,p(f^^^.) = CiM < C. 

To prove the invertibility of V$(u) = I + Vw in M? \ Vt^f we notice the following 
estimate for the first derivative of the solution to Laplace equation [19]: 

|Vw| < ||w||(^i(iB,3\f^^^) < C2||u||ci,i-3/p(r^.) < C2M < C, 

Therefore if M is chosen such that 

^' ^ ,na.{Cl.C2} ».5) 

V$(u) is an invertible matrix in M^. □ □ 

Remark 3.4. It follows from Lemma (3.3) that the matrix F(u) in Eq.(1.5) is well- 
defined, symmetric and positive definite. More precisely, we have that the maps F(u) (x) G 
anJ F(u)(a;) G C°^(]R'^ \ i^m/)- On the other hand, as a mapping from 
u G W'^'P{Qrnf) to F(u) G C°^(]R^ \ ^Imf), F(u) is infinitely differentiable with respect 
to u. In all what follows we will write F(u) and J{u) as F and J only, keeping in mind 
that they are u dependent. 
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4. Preliminary Results for the Poisson-Boltzmann Equation 

4. 1 . The Poisson-Boltzmann equation with Piola transformation. The rigorous anal- 
ysis and numerical approximation of solutions to the Poisson-Boltzmann equation (1.2) 
or its transformed version (1.4) are generally subject to three major difficulties: 1) the 
singular charge distribution, 2) the discontinuous dielectric constant on the molecular 
surface and 3) the strong exponential nonlinearities. However, it was recently demon- 
strated [2] that as far as the untransformed Poisson-Boltzmann equation (1.2) is con- 
cerned, some of these difficulties can be side-stepped by individually considering the 
singular and the regular components of the solution. Specifically, the potential solution 
is decomposed to be 

(j) = G + (f)'' = G + (f)^ + (f)'' (4.1) 



where the singular component 



(4.4) 



I 

is the solution of the Poisson equation 

N 

-V■{emVG)=pf■.= Y,(l^^ixi) in M^; (4.2) 

i 

while 0' is the linear component of the electrostatic potential which satisfies 

-V ■ (eV0') = - V ■ ((e - em)VG) in fi, 

, (4.3) 

(j) = g — G on (9fi, 

and the nonlinear component 0" solves 

-V ■ (eV0") + sinh(0" + 0' + G) = in 

0" = on dVL, 

where 

9 = y,<l~i 1 (4.5) 

is the boundary condition of the complete Poisson-Boltzmann equation (1.2). Such a 
decomposition scheme removes the point charge singularity from the original Poisson- 
Boltzmann and it was shown in [2] that the regular component of the electrostatic poten- 
tial 0*" = 0' + 0" belongs to (f2) although the entire solution G + 0'' does not. The 
most prominent advantage of this decomposition lies in the fact that the regular compo- 
nent represents the reaction potential field of the system, which can be directly used to 
compute the solvation energy and other associated important properties of the system. 
It is not necessary to solve the Poisson-Boltzmann equation twice, once with uniform 
vacuum dielectric constant and vanishing ionic strength and the other with real physical 
conditions, to obtain the reaction field [29]. As to be shown later on, the identification of 
this regular potential component as the reaction field also facilitates the analysis and the 
computation of the electrostatic forces. 

Applying the similar decomposition to the transformed Poisson-Boltzmann equation 
we get an equation for the singular component G: 

- V ■ (e^FVG) = Jpf in M^ (4.6) 



(4.8) 
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and an equation for the regular component (j)^': 

-V ■ (eFV0") + Jk^ smh(0" + G) = V ■ ((e - e^)FVG) in fi, 

(jf = g-G on dn. ^'^■'^^ 

We shall prove the existence of (jf in Eq.(4.7) and give its L°° bounds by individually 
considering the equation for the linear component 0': 

-V ■ (eFV0') = V ■ ((e - e^)FVG) in fi, 

4)" = g — G on dVt, 

and the equation for the nonlinear component 0": 

-V ■ (eFV0") + Jk^ sinh(0" + 0^ + G) = in fi, 

0" = O on ^"^'^^ 

As mentioned above, the functions G , cf)^ , and k in Eqs.(4.6) through (4.9) shall 
be interpreted as the compositions of the corresponding entries of these functions in 
untransformed equations (4.2) through (4.4) with the Piola transformation i.e., g = 

^($(x)), G = G($(x)), 0' = 0'($(x)), 0" = 0"(<l>(x)), = p-^($(x)), K = k(<I>(x)). 

4.2. Regularity and estimates for the singular solution component G. We first study 
the Eq. (4.6) for the singular component of electrostatic potential. We remark that the 
linear and nonlinear PB equations have the same singular component of the electrostatic 
potential. The solution of this singular component is the Green's function for the elliptic 
operator L defined by 

= -V ■ (e^FVw). (4.10) 

We shall use the following theorem [11] concerning the regularity and the estimate of the 
Green's function: 



Theorem 4.1. Let Q be an open set in M?. Suppose the elliptic operator 

d du . 



is uniformly elliptic and hounded, while the coefficients a^j satisfying 

\aij{x) - aij{y)\ < u{\x - y\) 
for any x,y E Q, and the non-decreasing function u{x) satisfies 

u{2t) < Kuj{t) for some K > Q and allt > 0, 

u{t) 



t 



-dt < oo. 



Then for the corresponding Green's function G the following six inequalities are true for 
any x,y E Q: 

(a) G{x, y) < K\x — y\~^, 

(b) G{x,y) < K6{x)\x-y\~\ 

(b) G{x,y)<K6{x)6{y)\x-y\-\ 

(d) \VMx,y)\<K\x~y\-\ 

(e) \VyG{x,y)\<K6{y)\x-y\-\ 

(f) \V,VyG{x,y)\<K\x-y\-^ 

where 6{y) = dist(?/, dQ) and the general constant K = K{aij, u, Q). 
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From this theorem we can derive the regularity of the Green's function of the operator 
(4.10). Indeed, by Sobolev embedding e^F G C°'^~^/p(]R^), therefore it satisfies the 
conditions on aij in this theorem provided that uj(t) = Kt^/^. We then conclude that the 
singular component of the electrostatic potential G E W^'°°{^1 \ Br{xi)). On the other 
hand, from Eq. (4.6) we know that G'($(u)(x))/J(xj) itself is the Green's function of 
operator (4.10) if F is generated by the Piola transformation according to (1.5) and J 
is the corresponding Jacobian. Thus the Green's function of differential operator (4.10) 
belongs to W'^''P{ri\Br{xi)) since it is the composition of the Green's function of Laplace 
operator, which is of C°°(i7\i?r(a;i)), and the Piola transformation, which is of W'^'^^^l). 
Higher regularity of G in can be derived thanks to the harmonic extension of u to 
\ ^Imf- In particular, because all charges are located in fi^/ and ^Imr the Poisson 
equation (4.6) appears a Laplace equation 

V(eFVG') = in fi„ 

hence G{x) G C°°(f2s), since Vlg is a smooth open domain and F G C°°(fis). 

In addition to the regularity of the Green's function, we have following estimates of G 
with respect to F and J. 

Lemma 4.2. For any given molecule the Green' s function G of operator (4.10) has esti- 
mates 

(a) ||G||^oo(n^) < C|| J||loo(j^). 

(b) l|VG|L^(n^)<C||J|U^(f,). 

If in addtion ||F — I||H'i.p(n) < Cf, \\J — l||vKi.p(n) < Cj for some constant Gf and Gj, 
then 

(c) IIGIlLPcan) < C'IIG'IIj^oo(q^). 

(d) ll^f o $||M/2-i/p,p(aQ) < Gg\\g\\w2,P{n^y 

(e) ll^f o $ - G\\]^2-i/p,p(^Qfi) < Gg\\g\\w2,p(n^) + GcWGHLocf^Q^y 

(f) ||FVG'||p^i-i/p,p(r) < Gr\\G\\ LOO (^n'^) for some set Q'^. 

Proof This ||^||j^oo(n^) is well defined since J is uniformly continuous in Us- To prove 
(a) and (b) we define qmax = niax{|gj|} and 

K K 

II VxG'j(a;, Xj) ll^oo(Q^) = ||G'i(x,a;i)||^oo(f^^,) = — 

where 5 is the smallest distance between x G dVt and singular charges at Xi. This smallest 
distance is related to the radii of atoms used in defining the molecular surface. In the 
sense of Connolly's molecular surface, 5 is simply the smallest van der Waals radius of 
the atoms which have contact surface [33]. We can therefore bound G and its gradient 
with 

ll^llL°°(n^) = II JQiGi\\Lao{fi^) < Nqmax\\J\\L°<^{fi^)\\Gi\\i^oc(Q^^-) 

i 

max 



6 

\'^G\\L°°(n^) = II 'Jli'^Gi\\L°°(ns) — ^^rnax\\J\\L°o(ns)\\'^xGi\\i^oo(Q^) 



(4.11) 



^(il)NKqmax 



52 

where is the total number of singular charges and || ^||Lcx)(nj,) is the maximum Jacobian 
onr. 
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The Statement (c) holds because dVl is also a piece of boundary of as shown in 
Fig.(l). To verify the statement (d), we noted that (7 o $ is the composition of g in 
Eq.(4.5), which is smooth in ^Is, and the mapping G W'^'^{^ls), i.e., 



^0$ = 



^-K\<;>{x)-<i>{xi)\ 



Following the estimate of the composite function in Sobolev space [12], we have the 
inequality 

\\g o ^\\w^-^/p'P{dn) < IkllvK^.pcfi,) 

< C(l + ||$||Loo(f^^))(l + ||$||vi/2,p(f^^))||5(||^2,p(n^) 

:= Cg\\g\\w2.p(ns) (4.13) 

with a constant Cg depending upon Here we choose to bound H^f o $||^y2-i/p,p(gj^) 

by \\g o $||^2,p(Q^) instead of \\g o $||vi/2,p(q) since the latter is not well defined due to the 
singular nature of g. 

The validity of inequalities (e) and (/) follows from the estimate of ||G||^2,p(j^'). This 

il'g is chosen such that ilg CC il'^. For example, we can choose to be the union of ilg, 
r, dVl, the domain 

= {x\x G Qrnf, dist{x,r) < -}, 

and the domain 

n+ = {x\x ^ n,dist(x,a(]) < -}. 
Applying the estimate to Eq.(4.6) in ^Ig we obtain 

||G'||p^2-i/p,p(aQ) < C\\G\\w2,p(ns) < C'(F) 11^11^^(5^^) < C(F)||G'||^oo(q'j 

:= (4.14) 

where the second inequality is a consequence of the estimate of the solution to — V • 
(eFVG) = in ^1'^. The coefficient Cq = C(F) depends on the ellipticity constants 
of F and its moduli of continuity on Vlg, hence is bounded as long as F is bounded. By 
combining Eqs.(4.14) and (4.13) we get (c). For the last estimate we notice 

||FVG||vi/i-i/p,p(r) < C|| [e]FVG||M/i,p(n^) 

< C\\F\\wi,PinJ^G\\w^.,^n.), (p > 3) 

< C'||F||vi^i,p(n^)||G||vy2,p(Q^) 

< C{F)\\F\\whP(^n^)\\G\\LP(n',) 
^ C'||F||vi^i,p(n^)||G'||L°°{ng 

:= Cr||G|Uoo(f,,). □ (4.15) 

□ 

Remark 4.3. ||G||Lo=(r2^^) can also be estimated by Eq.(4.11) if 5 is replaced by 5/2 and 

\\J\\L'^{n) is replaced by || J||L°°(nufi+)- 
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4.3. Regularity and estimates for the regular linearized solution component 0^ We 

consider an elliptic interface problem modified from the Poisson-Boltzmann equation 

-V ■ (eFVf) + Jn^f + G) = V ■ ((e - e„)FVG') + / in 

[0T = 0s-C = O, onr (4.16) 

(j)^' = g — G on dil, 

where [e] = — is the jump of dielectric constant and / G Lp{^1) is a given function. 
The equation for the regular potential solution of the linear Poisson-Boltzmann equation 
is a special case of (4.16) with / = 0. We remark that the regular component of the 
linear Poisson-Boltzmann equation in the absence of the Piola transformation represents 
a typical elliptic equation with discontinuous coefficients, for which Theorem (2.2) can 
be directly applied to get the existence and the estimate. In fact, the potential solution in 
this case is smooth in every subdomain (Proposition 1.4, [16]). When the Piola transfor- 
mation is incorporated, the coefficients of the Eq. (4.16) are not smooth and we have to 
rebuild the regularity and the estimate of the regular potential solution 0*", as summarized 
in the following theorem 

Theorem 4.4. There exists a unique solution cff of (4.16) in H^{Q). Moreover, there 
exists a positive constant C f such thatif\\F — V\\yYi,p^p^-) < Cj then (p'' belongs to >V^'^(f2) 
and the following estimate holds true 

WWw^'Pin) < G2 {\\G\\LP{n,) + \\f\\LP(n) + \\g - G\\w2-i/p,pi^gQ-)+ 

IIFVGII (4.17) 

Before proving this >V^'P estimate, we first establish a lemma concerning the esti- 
mate of a linear elliptic interface problem. 

Lemma 4.5. Let (f)^' solve 

-V-(eV0O + &0' = fir^^ 
[0^ = 0, onT 

M = 0, onr 

(f)^ = g on dQ, 

where e is a piecewise constant as defined for problem (4.16) and b > is a given real 
number, f{x) G LP{Q),ge H^iQ),p > 3. Then 

(4.18) 

Proof. The existence of unique solution 0'' G >V^'P c (i7) can be directly deduced 
from Theorem (2.2). We follow [2] and let cjf = cp^ + cj/' where 0' solves 

-V • (eV0^) = fmVL 



OonT, 
[e0y = Oonr, 
g on dVl, 



and 0" solves 



-V ■ (eV0") + 6(0" + 0') = Oinfi, 

[0"] = Oonr, 

[e4>l] = Oonr, 

0" = {)ondVt. 
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It is well known [20, 21] that 

II All 



\L° 



< C {\\f\\LP{n) + \\g\\m/2(^QQ)) 



while for 0" we claim that — ||0'||loo < ||0"||ioo < ||0'||ioo. To prove this assertion 
we define (f)t = max(0" — a, 0) where a = ||0'||loo. Then the trace Tr(0f) = hence 
0t G Hq by definition. Consider the weak formulation of the problem for 0" with 
test function 0t 

(eV0",V0i) + 6(0" + 0^0t) =0. 
Since 0* > wherever 0" > a, we have 

&(0" + 0', (j)t) = [ + (p')Mx + [ + <P')Mx > 0, 

and 

> (eV0^ V0t) = (eV(0" - «), V0i) = (eV0i, V0t) > 0. 

Thus V0t = 0, and 0t = or 0" < a in follows from the Poincare inequality. By 
defining (pt = min(0" + a, 0) and following the same procedure we can verify that 
0" > —a. The lemma shall be finally proved by combining the estimates of 0' and 

0". □ □ 

Proof, of Theorem (4.4). Consider the general weak formulation of the elliptic equation 
in problem (4.16), i.e., find 0^ = m e H^{Q) such that A{u,v) = F{v),\fv G H^{n) 
where 



A{u, v) = (FVuVf + JK^uv)dx, 
Jn 

F{v) = [ {V ■iie-6m)FVG)-JK^G + f)dx-A{g-G,v). 
Jn 

We shall apply the Lax-Milgram theorem to obtain the existence and the uniqueness of 
a weak solution 0*^ G H^{ri) to (4.16). Hence we must show that F(-) is bounded, and 
y4(-, •) is bounded and coercive with the assumptions on the coefficient matrix F and the 
Jacobian J. Consider the bilinear form A(-, ■). The Piola transform matrix F is positive 
definite, hence FVf ■ Vf > 7I Vf p for some 7 > 0. This inequality and the positiveness 
of Jacobian J give 

A{v,v) = (FVv -Vv + JK'^v'^)dx > {-flVvl"^ + jK^v^)dx > X\u\H^n) 
Jq Jn 

= 7 Qbl^MQ) + ll^lmin)^ > 7 (^-^IHhn) + l\^\m(n)^ 

> m(^\\v\\l2^^) + \v\l,(^^)^ =m\\v\\li(^^^, (4.19) 



where in the second inequality we applied the Poincare inequality with constant 6. Thus 
we verified that A(-, ■) is coercive, with coercivity constant m = min{7/(26'^), 7/2}. 
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On the other hand, 

\A{u,v)\ = \ {FVu-Vv\ + JK\v)dx\ (4.20) 
Jn 

— / |Fjj-Dj'uDjt'|(ix + / \jK'^uv\dx 

< ^ \\Fij\\Lo.^n)\\DiuDjv\\Li{fi) + fi^\\J\\L°°\\uv\\Li(^n) 

id 

— l|Fij||L°°||M||Hi(n)||'w||Hi(f7) + /t^||-'^||L°°||M||L2(n)||'i^||L2n 
id 

< Ki\\u\\Hi{n)\\v\\m{n) (4.21) 
which proves that A{-,-) is bounded with constant Ki = ||Fjj||i^oo + J||/^cx). This 

id 

constant Ki is finite because F, J belong to W^'P{il) which is compactly embedded in 
C°(n) forp > 3. 

In order to apply the Lax-Milgram theorem it remains to show that F(v) is bounded 
onH^{n). We have 



\F(v)\ < 



[ {e-em)FVGvdx\+ [ \Jk'^Gv + fv\dx + \A{g - G,v) 
Jq Jq 

)FVGvdx+ / {es - e„,)FVGvdx\ + 

J S^TTT, J fig 



Jk'^Gv + fv\dx + \A{g - G,v)\ 
= I / {ts - em)FVGvdx\ + / \Jk^Gv + fv\dx + \A{g - G,v)\ 

< [ lies - em)FVGv\dx + [ \Jk'^Gv + fv\dx + \A{g - G,v)\ 

< [e]\\FVG\\L2^n)\\v\\LHn) + in^JG\\L2^n) + \\f\\LHn))\\v\\LHn) + 

= ([e]||FVG|U.(f,) + K'\\JG\\L2^n) + ll/IU^cn) 

- G\\H^(n)) \\v\\H^n) 
= K2\\v\\Hi(n), 

hence F(-) is a bounded linear functional on (i7). 

We now proceed to show the regularity result and the estimate of (f)^ following the sim- 
ilar iterative technique in [6]. For this purpose we introduce a sequence {0^} generated 
by 

- V ■ (eV0^) + JkV^v = V ■{{e-em)FVG)-JK^G (4.22) 

+V-(e(F-I)V0^_i)infi, 
[0^] = Oonr, (4.23) 
0^ = g — G on dil, 
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and prove that 0^ G >V^'^(r2) and 0^ converges to the unique solution (jf of (4.16) in 
W^'^(f2) as — 7- oo. The first term 0q of the sequence solves 

- V ■ (eV0S) + Jfi:205 = V ■{{e-e^)¥VG)- JK^GmVt 

Wo] = 0, onr (4.24) 
0Q = g — G on dVl, 

therefore it belongs to >V^'^(i7) according to Theorem (2.2). Suppose now that G 
W'^'P{n), then V0^_i G W1'P((^) and V ■ (e(F - I)V0^_i) G LP(^]) following from 
Lemma (2.1). Thus problem (4.23) also has a unique solution 0^ G >V^'^(^7) for all 
integer according to Theorem (2.2). To prove that converges to the unique solution 
0*" of problem (4.16), we estimate ||0^ — 0^_i||w2,p(f^) and show it is decreasing as 
— 7- oo. By subtracting the equations in (4.16) for from those for A^ — 1 we obtain a 
problem for — (j)'^-!- Applying Theorem (2.2) again we know that this problem has 
a unique solution in W'^'^^^l) which has an estimate 

Un - fN^ilU^nn) < C(||V-(6(F-I)V(0^_i-0^^_2))IU.(o) 

+ Un - <pN-l\\LP{n)) 

< C||V ■ (6(F - I)V(0^_i - 0^,_2))IU.(C), 

< C||F - I||vyi,P(f7)||0^_i - 0^_2llw2.p(n), (4.25) 

where in the second inequality we applied Lemma (4.5) to the problem for (0^ — 0^„i), 
and the generic constant C is independent of A^, F. Therefore if the constant C/ in the 
assumption of the theorem is chosen such that GGf = k < 1 then ||0^ — 0^_^||yv;2,p(Q) 
is decreasing with respect to A^ hence the sequence converges to a unique element (p^ 
in >V^'^(fi). Letting A^ — > oo we can observe that (j)^ is the unique solution of problem 
(4.16), meaning 0^' = 0''. 

The estimate of (p^' is obtained by estimating 0^ and passing A^ to oo. We notice that 
07V = 07V - 0iv-i + 07V-1 - 0iv-2 + • ■ ■ + 0S' hence 

1 _ k^-^ 

||0^||w2.P < _ ^ 1101 - 0Sllw2,P + ||0o||h;2,p 

< C2(||G'||LP(ns) + \\g — G'||,4/2-i/p,p(f^) 
+ ||FVG||v7i-i/p,p(r) as A^ ^ oo, 

where both Theorem (2.2) and Lemma (4.5) are applied to the problem of (pl — 0q and 
the problem of 0q to get the desired bounds with respect to the W'^'^ and norms, and 
C2 absorbs k and all the generic constants involved in these bounds. □ □ 

4.4. Regularity and estimates for the regular nonlinear solution component 0''. For 
the nonlinear Poisson-Boltzmann equation, the regular component 0*" of its potential 
solution solves 

-V ■ (eFV0") + Jk^ sinh(0" + G) = V ■ ((e - e,n)FVG) in 

[0T=0^-0:„ = Oonr (4.26) 
0^ = ^ - G on dn. 

The appearance of the nonlinear function sinh(a;) complicates the establishment of the 
existence of (p^. In particular, the Lax-Milgram Theorem is not applicable to problem 
(4.26). Instead we define a energy functional based on the weak formulation of (4.26) 
and show that the unique minimizer of this energy functional is the unique solution of 
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(4.26). On the other hand, the establishment of the regularity and W^'^ estimate of cjf' for 
(4.26) is simplified thanks to Theorem (4.4). 
We start with the weak formulation of (4.26): 

Find (jf eM = {ve {VL)\e\ e-'' G ^^(fi), and v = ^ - G on dVl}, 

such that 

A{(t>\v) + {B{(^n,v) + {fc^v) = 0, Vi; G Hl{n), (4.27) 

where 

A{<f, v) = (eFV0^ Vv), {B{f), v) = {Jk^ sinh(0'^ + G),v), 

ifcv) = [ {e-e^)FVG-Vv. 
Jn 

We also use fc to denote the function [e]F(u) VG ■ n on the dielectric boundary T, since 

ifcv) = {[e]FVG-n,v), (4.28) 

where [e] = — is the jump in e on F. Based on this weak formulation we define an 
energy on M: 

E{w) = I -FVw ■ Vw + Jk^ cosh(w; + G) + {fc, w). (4.29) 
Jn 2 

The weak solution of Eq.(4.7) can be characterized as the minimizer of this energy func- 
tional. This equivalence and the existence of this minimizer are due to the following four 
simple lemmas. For the proof of these lemmas we refer to [2]. 

Lemma 4.6. Ifu is the solution of the optimization problem, i.e., 

E{u) = inf E{w), 

then u is the solution of ( 4. 7). 

Lemma 4.7. Let F[u) he afunctional defined on M, if 

(1) M is weakly sequential compact, and 

(2) F is weakly lower semi-continuous on M, 

then there exists u & M such that 

Fiu) = inf F(w). 

w&M 

Lemma 4.8. The following results hold true 

(1) Let V be a reflective Banach space. The set M := {v E V\\\v\\ < tq} is weakly 
sequential compact. 

(2) if lim F{v) = oo, then inf F{w) = inf F{w) 

Lemma 4.9. If F is a convex functional on a convex set M and F is Gateaux differen- 
tiable, then F is w.l.s.c. on M. 

The existence and the uniqueness of the weak solution to (4.7) can be established 
using these lemmas. The following lemma establishes the existence of the minimizer of 
the energy E{w). 

Theorem 4.10. There exists a unique u E M C. H^{Q) such that 

Eiu) = inf E(w). 

w&M 
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Proof. The differentiability of E{w) follows its definition. Actually we have 

{DE{u),v) = A{u,v) + {B{u),v) + {fc^v). 

The minimizer of E{w) exists if we can prove that 

(1) M is a convex set 

(2) E is convex on M 

(3) lim E{w) = oo 

It is easy to verify (1). The convexity of FVw ■ Vw follows from the fact that 

< FV(7w) ■ V(7w) = 7^FVw ■ Vw < --fFVw ■ Vw 

for any < 7 < 1 since F is positive definite. The convexity of cosh(iy + G) follows 
from the convexity of cosh(a;) directly. Actually E{w) is strictly convex. To prove (3) 
we only need to show that 

Eiw) > Cie,K,F)\\w\\l,^^^+CiG,g). (4.30) 

We notice that cosh(a;) > 1 and 

<fG,W> < e^||FVG||L2(f^^)||Vu7||L2(f^^) 

< es\\F\\L2^a,)\\VG\\L2(n,)\\Vw\\L2(ns) 

< ^(l|VG||i.(^^) + ||Vu;||i.(^^)), 

where the matrix norm 



7 = l|F||L2(n,) = sup \\Fv\\L2(n^), 
feM,||i)||^2(n)=i 



is finite because F is continuous. Therefore 



E{w) > - [ eFVw - Vw - \ < few > \ 
2 Jn 

= ll e„|Vt.p-^||VG|| 



> C(6,7)||V^||i.(^)-^||VGL.(^^). 



!^IIV7^l|2 

2 

The inequality (4.30) follows from the equivalence of || Vw||L2(f^) and on set 

M. The uniqueness of the minimizer of E{w) comes from the strict convexity of E. □ 

□ 

Theorem 4.11. There exists a unique solution cff of (4.16) in H^{Q). Moreover, there 
exist constants Gi, G2 and C3 such that cj)"^ is bounded by 

<Gi + C2II J||L°°(f7) + C3||-'^||L°°(n)||F||tyi,p(n^). (4.31) 

Proof. The existence of the solution (p''' in H^(n) has been proved by theorem (4.10) 
and its four lemmas. It remains to verify the L°° bounds of 0''. Let 0^ = 0' + 0" be 
decomposed into a linear component 0^ and a nonlinear component. The linear com- 
ponent 0' satisfies Eq.(4.8). The existence of a weak solution 0' G Hq^VI) follows that 
V ■ ((e — em)FVG') is an operator in H^^(Q) [8]. It is well known that in general 

Ci||0'|U^ < U'Wm < C2i\\9 - G||HV2(ac) + WfaWm/^r))- (4-32) 
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To estimate the nonlinear component we follow [2] and define: 

ol = argmax{ Jk^ sinh(c + sup 0' + sup G) < 0}, 

(3' = argmin{ Jk^ sinh(c + inf 0^+ inf G) > 0}, 

a = min(a',0), 
(] = max(/3',0). 

It follows from the monotonicity of sinh(x) that 

P = ||0'IU-(f^,) + ||G|Uoo(o,), (4.33) 
a = -13. (4.34) 

We will show that a and (3 are the lower and upper bounds of the nonlinear compo- 
nent 0" of the weak solution to (4.26), following the similar procedure as that used in 
proving Lemma (4.5). 
Define 

(j)t = max(0" -13,0), 

then Tr0f = since 0" G Hq(^1) and /3 > by definition. Therefore (pt E Hq{^1) and 
satisfies the weak formulation of Eq.(4.9): 

(eFV0", V(j)t) + {Jk^ sinh(0" + 0' + G),(j)t) = 0. 

Since (f)t > wherever 0" > j3, we have 

J{u)k'^ sinh(0" + 0' + G) > J(u)k2 sinh(/3 + inf 0' + inf G) > 0. 

x&Qg x(^Qs 

Therefore 

> (eFV0", V0i) = (eFV(0" - V0t) = eFV0t • V0t > 0, 

where the last inequality holds true since F is positive definite. Hence V0t = 0, and 
0t = or 0" < /3 in i7 follows from the Poincare inequality. This establishes the upper 
bound of 0". By changing (pt to be min(0" + a, 0) we can also prove that a is the lower 
bound. 

Combining the estimates for 0' and 0" we finally obtain the L°° estimate of the regular 
component 0^: 

W + rho^m < ||0^1|Loo(n) + ||0l|L-(Q) 

< ||0'||l°°(!^) + ||0'||L°°(n,) + 11^111,00(5^^) 

< 2{Gg + Ccll-'^llLa^cn) + C'/Gl|F||iyi.P(n3)l|-'^IU°°(!^) 

max 

6 

= Ci + C2II J||L°=(f7) + C'3||-^||L°°(n)||F||vi/i.p(ns)- n 



□ 



We are now able to examine the regularity results and the estimate of 0*^ in W^'^. 



Theorem 4.12. If\\F — IllH'i.p(n) < Cf then the unique solution (jf of (4.26) belongs to 
W'^'^iVL) and the following estimate holds 

\W\\w^'P(a) < G {\\G\\Lv{ns) + \\g — G'||vi/2-i/p,P(aQ) + ||FVG||^yi-i/p,p(r)) . (4.35) 
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Proof. It is noticed that the problem (4.26) can be written as a form similar to its linear 
counterpart (4.16) 

-V ■{eF\7(lf) + JK\(jf + G) = V- ((e-e^)FVG) 

-jK\smh{(j)'- + G)- (0" + G))mQ 
[f] = 0, on r 
(J)"- = g-Gondfl. 

According to Theorem (4.4), the W^'^ regularity of (p^ directly follows from the facts that 
V-((e— em)FVG) represents an interface condition in W'^^^^^''p{^1) and that J/t^(sinh(0''+ 
G) — {(jf + (?))€ In the mean time, we have an estimate 

110' ||w2'P(n) < C {\\G\\lp{q.s) + Wd — G^||iy2-i/p,p(an) + l|FVG||vyi-i/p,P(r)) • □ 

□ 

5. An Electrostatic Force Model and Some Estimates 

For the untransformed nonlinear Poisson-Boltzmann equation (1.2) the electrostatic 
energy of the system is defined [27, 13, 14] to be 

E = I W<P - \eiy<Pf - «:'(cosh(0) - l)x]dx, (5.1) 

where the characteristic function x = 1 in and is in molecules Vtmf, ^mr- This 
energy is very similar to the energy functional defined in Eq.(4.29), and any potential 
function minimizing (4.29) is also the minimizer of this electrostatic energy because 
cosh(a;) > 1. The function cosh(0) — 1 describes the physical fact that the total electro- 
static energy is zero when is everywhere zero. The three terms in this energy represent 
three types of energy densities, namely, the Coulomb energy, the electrostatic stress en- 
ergy and the osmotic stress energy of the mobile ions. Based on this energy function, 
the following density function of the force exerted on the molecule was derived [14] by 
using a variational derivation method: 

f = p-^E - ^|EpVe - /t2(cosh(0) - l)Vx, (5.2) 

where the three terms correspond to the Coulomb force, dielectric pressure and the ionic 
pressure, respectively. The last two boundary forces are always in the normal direction 
of the molecular surface because of the gradients of e and the characteristic function x- 
The electric force defined in (5.2) is physically justifiable, and can be converted into a 
form identical to the Maxwell stress tensor(MST) [13, 32]. The MST describes the vol- 
ume force density in a linear dielectric, and has been widely utilized in dielectrophoretic 
force and electrorotational torque calculations of colloids, macromolecules and biologi- 
cal cells in continuous external electric field [32]. In the context of interactions between 
singular charges distribution and resulting singular electric field, refinements are neces- 
sary to make this force model computationally more tractable. Below we will discuss the 
treatments of its three components. 

The first term in Eq. (5.2) might appear misleading because of the multiplication of 
two singular functions, and E, in its expression. We therefore would emphasis that 
at a singular change Xi the electric potential field multiplied with in Eq.(5.1) shall be 
interpreted as the summation of reaction potential field cj/ , i.e., the regular component of 
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the potential solution, and the Coulomb potential induced by all other singular charges 
[14]: 

p^E = 5^g,5(x,)E := V U^'(a;) + ^.(x) J 

= J^g.V U^(x,) + ^G,(x,) J (5.3) 

This verifies that the force exerted at each charged atom is finite. The eliminated term 
Gi{xi) corresponds to the self-energy of the singular charges [14]. 

Nevertheless, the body force density p-^E itself is still unbounded at the center of every 
charged atom where the charge density is singular, indicating that this body force density 
does not belongs to hence does not fit the assumption on the body force in 

Theorem (3.1). An alternative model is therefore necessary to regularize these singular 
body forces to ensure the solvability of the elasticity equation. In this study, the singular 
the body force density is modeled by a Gaussian function 

fb = 5^a,e-(^-^»)'/'^»n„ (5.4) 

i 

where the unit normal vector is aligned with the corresponding gradient in Eq.(5.3); the 
decay parameter a-i is chosen such that 

^aie-^?/"^ = 5 (5.5) 

i 

for a given sufficiently small number 5, and Ri is the van der Waals' radius of atom i. 
This means that the Gaussian function is essentially compact supported in its associated 
atom. The prefactor ai is determined by the conservation of force in each atom: 

I aie-(^-^')'/"'da; = q^[<if{xi) + ^ G^{x,)] = iTiai j " r^e-^'/^'dr. (5.6) 

J atonii J 

The body force fb modeled by this Gaussian is uniformly continuous in Vtmf and belongs 
to L^iytmf) for any p > 0. Moreover, the lemma below proves that the difference of two 
continuous body force densities also belongs to L'^iVtrnf), and is small if the difference 
between two total body forces which they approximate is small. 

Lemma 5.1. Let Ai,A2 be two given numbers and \Ai\ < P,\A2\ < P for some P. Let 

f. = aj.e-(^-^o)'/'^^ such that j a^e-^^-^'^^'/'^^c/x = Aj for j = 1,2, 

J x—xo<R 

where aj, aj are determined from Eqs.(5.6,5.5) for the same atom centered at xq and of 
radius R. Then if\Ai — A2\< 5' for some 5' > 0, we have 

\fi-f2\'dx<C5' (5.7) 

x—xo<R 

for some constant C depending only on R and P. 

Proof The prefactor a and the decay rate cr are uniformly continuous functions of A for 

\A\ < Pif 
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is the approximation of A as defined by the lemma. But then there exists a constant C 
depending on the derivatives of a and a with respect to A such that | fi (x) — /2 (a;) | < C5' 
if |Ai — < 5'. The conclusion of the lemma follows directly. □ □ 

The last two terms in Eq.(5.2) represent the electrostatic surface forces on the mole- 
cule. It is worth noting that the second term is not well defined and is computationally 
intractable if there is no dielectric boundary smoothing, due to the discontinuous electric 
field E on the molecular surface indicated by the interface condition 

emV0m ■ n = e^V^s ■ n or e^E^ ■ n = e^E^ ■ n. 

To remove this ambiguity we consider a infinitesimal displacement h of the molecular 
surface in its out normal direction, see Fig. (2). The change of the electrostatic stress 
energy due to this small displacement is the work done by the dielectric pressure along 
this displacement: 

^e|Ep(ix- / --e\E\'^dx = [ --(es|Es|^ - e„|E„p)(ix 



feds. 



This suggests the dielectric force density fg is essentially the difference between — |es|E<j 
IE I 

2 "-rra I \ 



and — I Cm I Em P on the dielectric interface, i.e.. 



fe = -^(e,|E,|2-e„|E„nn. (5.8) 

By combining definitions (5.2), (5.4) and (5.8) we would obtain a complete model of the 
electrostatic body force and surface force: 

-ix-x,f/a. 



(5.9) 



fs = -^(e,|E,|2-em|Emnn-/€2(cosh(0)-l)n. (5.10) 

Remark 5.2. In the sequel we will estimate the term ||fs|| vKi-i/p.p(r/)- Although the term 
||Es||^i-i/p,p(P^) can be directly related to \\4>s\\w^'P{vls) since the latter term is bounded 
in VLs, one can not estimate \\Em\\w'^-^/P'P{rf) similarly by relating it with \\(t)m\\w'^'P{n,y^f) 
because 0^ contains singularities and hence is unbounded in flmf- Instead we follow 
the procedure in the proof of (e,f) in Lemma (4.2) and eventually estimate this trace norm 
of E in which does not contain potential singularities; the details are omitted due to 
similarity of these two proofs. 

Remark 5.3. The surface force definition presented in Eq.(5.10) applies only to the dis- 
continuous dielectric model as adopted in this study. In the continuous dielectric models, 
which are also widely used for in the implicit solvent simulations, different surface force 
definition will be derived [10]. However, the analysis on the electrostatic forces given in 
the below is also applicable to general surface force function = /s(Es,Em,0), and 
might be simplified if electrical field E is continuous, i.e., e is continuous on T. 

The electrostatic forces defined in Eq.(5.9) and Eq.(5.10) are also subject to the Pi- 
ola transformation. Moreover these forces can not be directly supplied to the elasticity 
equation; only the forces relative to a reference state can be supplied. This is because 
a molecule is in an equilibrium state and has no elastic deformation if the electrostatic 
potential is induced only by the molecule itself and the solvent with physiological ionic 
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'S 



Figure 2. Displacement of the molecular surface Fj. The solid black 
line is the surface before displacement and the dashed red line is the sur- 
face after displacement. The new solvent region Q'^ is plus the strip 
between two surfaces; The strip subtracted from the equals the new 
solute region il'^. 

strength, in the absence of interactions with other molecules. We refer to this state as the 
free state and use it as the reference state. The net body force or the net surface force 
is therefore defined to be the difference between that for a molecule in non-free state 
and that for the same molecule in the free state. To abuse the notation these differences 
are still referred to as the body force and the surface force, and are denoted by f;, and 
respectively: 



where f^o and f^o are the body force and the surface force in the free state, and are constant 
vector fields for any given macromolecule. 

Physically, these two forces f;, and shall be vanishing if there is no change of ionic 
strength and no additional molecules present, and will be small for small change of ionic 
strength and weakly interacting additional molecules. To reflect this physical reality and 
to facilitate the mathematical analysis, we decompose (into four steps) the transition from 
the original single deformable molecule immersed in aqueous solvent with physiological 
ionic strength to the final system with added rigid molecules, varied ionic strength and 
deformed molecules. In the first step, we change only the solvent from physiological 
ionic strength to the target strength, and assume that the molecule ri„if does not have a 
conformational change although the net electrostatic force is not zero due to this change 
of ionic strength. The electrostatic potential and forces at the end of the first perturbation 
are denoted by 01 and f;,i, f^i, respectively. In the second step, we alter the dielectric 
constant in the smooth domain fi^r from to em- This low dielectric space represents 
the empty interior of the added molecules. The electrostatic potential and forces after 
the second step are respectively denoted by 02 and ft2, is2- In the third step we place 
the singular charges into and define the electrostatic potential and forces to be 03 
and fb3, fs3. In the last step we allow the Poisson-Boltzmann equation to couple with 
the elastic deformation so that the system will arrive at the final state with electrostatic 
potential and forces f;,, f^. We write the net body force f;, and the net surface force as 
the summation of their four components 



is 




(5.11) 
(5.12) 



ft 



(f, - U,) + (ffe3 - f62) + (f62 - fw) + (f61 - fbo) 
(f, - f,3) + (f.3 - f.2) + (fs2 - isl) + (fsl - f.o) 



(5.13) 
(5.14) 
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corresponding to the above decomposition, and estimate these components individually. 



5.1. The surface force due to changing ionic strength. The electrostatic potential 0o 
of the system in the free state is given by 

- V ■ (eV</)o) + Klsmh.{(j)o) = ^gi5(a; - Xj), (5.15) 

i 

while the electrostatic potential 0i after changing of the ionic strength satisfies 

- V ■ (eV0i) + ft^sinh(0i) = "^qiS^x - Xi). (5.16) 

i 

By subtracting Eq.(5.15) from Eq.(5.16) we get 

- V ■ (eV0) + (k^ - kI) cosh(O0 = (f^l - k^) sinh(0o) in n, (5.17) 

where = 0i — 0o and ^{x) G (min{0i(x), 0o(a:;)}, max{0i(a;), 0o(a;)}) is a function 
between 0i and 0o satisfying the Cauchy mean value theorem 

sinh(0i) = sinh(0o) + cosh(^)(0i - 0o). 

We note that the singular charges disappear in Eq.(5.17), and hence (p E H^{n) and 
is also in C°° in ^Imf and \ ^Imf- Moreover, following Theorem (2.2) we have the 
following W^'^ estimate for 0: 

||0||w2>f(n) < C (\\(f)\\LPin) + ||(/to - '^^) sinh(0o)||Lp(Q) + \\G\\Lp{dn)) , (5.18) 

where 

p~K\x-Xi\ _ „-Ko\x-Xi\ „-K\x-Xi\ 

G = y2- ^ ^-{n-Ko)y2- (5-19) 

I ' ' I 

is the boundary condition of on dVt, and is the difference of boundary values of 0i and 
00- The approximation in Eq. (5. 19) is well defined for small {k — kq). On the other hand. 
Lemma (4.5) says that ||0||LP(r2) itself can be estimate by 

||0||L.(n) < C||0|U-(n) < C - n^) sinh(0o)IU2(Q) + . (5.20) 

By combining Eqs. (5.18), (5.19) and (5.20) we get 

1101 - 0o||vK2.f(n) < C'I'to - (5.21) 

We now proceed to estimate the changes of electrostatic forces — f;,o, fsi — fso- The 
body force change 

llffoi - fM||LP(n„/) < C |0(xi)| < C|/to - k| (5.22) 
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follows from Lemma (5.1). On the other hand, 

fsi-fso = -^e^dV^i^P - iV^o^Hn + ^e„(|V0imp - iV^ornHn 

- (fi;^(cosh(0is) - 1) - Ko(cosh(0os) - 1)) n 
= -^es(V0is - V0os) ■ (V0is + V0os)n 

- /t^(cosh(0is) — cosh(0os))n 

- (k^ - Kg) cosh(0os)n 

We note that with the mean value theorem, ft^(cosh(0is) — cosh(0os)) can be related 
to the change of ionic strength as 

K^(cosh(0ij - cosh(0os)) = smh{C){(pis - 

Moreover, we can not bound the trace norm of | V0imP — | V0omP by its Sobolev norm 
in subdomain ^Imf where the singularities of the potential are located. Instead we follow 
the remark of Eq. (5.10) and estimate this term in subdomain fij. Thus the surface force 
change in the first perturbation step can be estimated as 



|fsi — iso\\w^-^/P:P{rf) < ~ 0o||iy2,p(j^^)||0i + 0o||vy2>p(na) 

+ 1101 - 0o|lvy2,p(Q-)||0i + 0o|lvF2.p(nj) 
+ sinh(^')lll0i - 0o|ki.p(n.) 



+ K 



i^o\Hos\\w^^p(ns) + \^'^ - ^ol) 

< C\k - Kol (^1101 + (po\\w^-p{ns) + 1101 + 0o||iyi,P(Qj) 

+ K^l sinh(^')l + + '«o)(||0olki>f{Q,) + 1) 

< C,(k)|k-ko|, (5.23) 

where Lemma (2.1) is applied to estimate the norm of the products of two W^''^ functions 

V01 - V0O and V0i + V0o. 

5.2. The surface force due to adding a low dielectric constant cavity. Although the 
variation of ionic strength will change the electrostatic potential of the system, the mag- 
nitude of potential change is usually smaller than that induced by adding molecules to 
the system. By adding molecules to the system we will not only have the additional sin- 
gular charges but also expand a cavity of low dielectric constant in the solvent. These 
two effects will be considered separately, and this subsection estimates only the change 
of potential and forces due to the additional cavity of low dielectric constant. The effect 
of added charges will be analyzed in the next subsection. 

The electrostatic potential 02 with an additional low dielectric cavity in the domain is 
described by 

- V ■ (eV02) + K^ sinh(02) = ^ qi5{x - Xi) (5.24) 
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Figure 3. Illustration of adding rigid molecule(s). Left: Before adding 
rigid molecules the domain fi^r is occupied by solvent and hence has 
dielectric constant e^. Right: After adding molecules the domain ilmr has 
low dielectric constant tm- 



with the same boundary conditions as Eq.(5.16). Here the dielectric constant e and ionic 
strength k are different from those in Eq.(5.16), and thus the subtraction of Eq.(5.16) 
from Eq.(5.24) shall be individually conducted in fim/, ^mr and fi^ to give the following 
three equations: 

-V ■ (e„V(02 - 0i)) = in n^j, 
-V ■ (e,V((/)2 - + K^(sinh((/)2) - sinh(0i) = in 1]„ 
-V ■ (emV02) + V(esV0i) - sinh(0i) = in Qrnr- 

By assembling these three equations we get a complete equation for = 02 — 0i in ^l: 

- V ■ (eV0) + K2(sinh(02) - sinh(0i)) = — smh{(j)i) , (5.25) 

where e is the same as that in Eq.(5.24), and the right-hand side is vanishing in and 
^Imf- This function (non- vanishing only in ilmr) is equivalent to —V ■ ((e^ — em)V0i) + 
ft^ sinh(0i) since 

— V ■ (esV0i) + sinh(0i) = in 

As before we notice that K^(sinh(02) — sinh(0i)) can be related to cosh(^)0 with a 
smooth function ^ bounded by (pi and 02, and therefore in Eq.(5.25) satisfies an esti- 
mate of the form 

U\\w2-p{n) < C |^||0||LP(n) + ||y^fi:2sinh(0i)||LP(n_)^ 

< C||sinh(0i)|U.(f,_) 

< C|| sinh(0i)||Loo(f^) ■ 14,^ 

which follows from Theorem (4.4) and Lemma (4.5), considering that Eq.(5.24) has a 
vanishing boundary condition. Here Vmi- is the volume of ^Imr, suggesting that ||0|| >v2 p(n) 
can be made arbitrarily small by reducing the volume of fimr- 

The change electrostatic body force induced by this additional low dielectric cavity 
can be estimated as 

< cY,\4>{x^)\<CVmr. (5.26) 
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The surface force change is 

fs2-fsl = -^es(|V02^r - IV^lsHn + ^em(|V02mP - iV^imHn 

— /t^(cosh(02s) — cosh(0is)n 
= -^e,(V02s- V01,) ■ (V02S + V0i,)n 

+ ^em(V02m - V0i™,) " (V02m + V0i„^)n 

-K^(cosh(02s) - cosh(0is))n, 
and thus can be estimated by 

||fs2 — fsl|lH/i-Vp,p(rj) ^ C'(^ll02 — 0l||l4'i>P(n,)||02 + 

+ 1102 - 0l|liyi.P(Q7)||02 + ^lllw-i^pCQD 

+ sinh(^')lll02 - 

< C||02 - 0l||w2>P(n) < C* ■ Knr, (5-27) 

following from the similar arguments in last subsection for estimating f^i — f^o- 

5.3. The surface force due to additional singular charges. In this subsection we will 
consider the change of electrostatic potential and force caused by singular charges placed 
in the low dielectric space Vtmr- The low dielectric space Vtmr with these charges com- 
pletely models the rigid molecule which is expected to interact with a flexible molecule 
VLmf- The electrostatic potential field after this third perturbation step satisfies the fol- 
lowing equation 

- V ■ (eV03) + sinh(03) = q,5{xi) + ^ qj5{xj), (5.28) 

« i 

while the change of potential, = 03 — 02 is the solution of the equation 

Nr 

- V ■ (eV0) + cosh(O0 = qjS{xj), (5.29) 

j 

which is obtained by subtracting Eq.(5.24) from Eq.(5.28). Here ^(x) is a smooth func- 
tion defined by the mean value expansion sinh(03) = sinh(02) + cosh(,^)(03 — 02). To 
facilitate the regularity analysis of we define its singular component G, which solves 

Nr 

-V-(e„VG) = ^g,(5(x,) (5.30) 

j 

and its regular component (p^, which is the solution of 

- V ■ (eV0'^) + cosh(O0'' = V ■ ((e - e^)VG') - cosh(OG'. (5.31) 

It shall be noted that V ■ ((e — em)VG) is nonzero only on the molecular surfaces Tf and 
r^, and can be represented as an interface condition (cg — em)VG' ■ n on each of these 
two molecular surfaces similar to that in Eq.(4.28). We notice that 

Nr 

^(^) = E -T^ 1 (5.32) 



I X Xj I 
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and is of C°° wherever away from any of Xj, hence of C°°(fis), and thus so is — (e^ — 
ej„)V(5 ■ n on and F^. The W^'^ estimate of in Vlg says that 

w2.p(n) < (^/t^ll cosh(^)G||LP(Q^) + IKe^ - e„)VG||t^i-i/p,p(r^) 

+ 11 (cs - em)VG'||,yi-i/p,p(p^) + ||5'llw'i-i/p,p(af7) + ||0lUp(fi)j 
< C (^\\G\\lp + IIGIlH/a.Pcn,) + \\9\\w^>p(n,)) ^ as ^ (5.33) 



where 



9 



Eqj — I r on dVL 



3 



is the boundary condition of Eq.(5.29), and is the difference of boundary conditions of 
Eq.(5.24)andEq.(5.28). 

We now analyze the change of the electrostatic forces due to the inclusion of additional 
singular charges. For body force we have 

\%z-%2\lv(si^,;) < C^|0(xi) + ^Gj(xi)| -4 0asgj- ^0, (5.34) 

and for the surface force we know 

- Je,(|V03.r - |V(^2.r)n + Je„(|V<^3m|' - |V02mr)n 



Is3 ~ Is2 — 

— (k^ cosh(03s) — cosh(02s)n 
= -^e,(V03, - V02s) ■ + V02s)n 

+ ^e„(V03m - V02m) " C^4>3m + V02m)n 

-K^(cosh(03s) - cosh(02s))n. 
This surface force difference can then be estimated by 



||fs3 — fs2||wi-i/p,p(r/) ^ C\^\\(p3s — (P2s\\W^-P{Qs)\\'f3s + (P2s\\W2.P{n^) 

+ \\(p3m — 4>2m\\w'2'P{n-)\\4>3m + 4>2m\\w'^^P{n7) 

sinh(^')lll03s - 02^lk2.p(n,)) 

< C'^ll'/'s - 02||w2.p{n,) + 1103 - 02||n'2>p{no) 

as qj (5.35) 
where the constant C depends on the W"^'^ norm of 02, 03, and therefore is bounded if 



^ 



P2,(P3 



are bounded. 



5.4. The surface force due to molecular conformational change. We now consider 
the change of electrostatic potential and surface forces induced by elastic displacement. 
By subtracting Eq.(5.28) from Eq.(1.4) and with a few algebraic manipulations we get 
the governing equation for = — 03: 

Nf+Nr 

- V- (eFV0) + J/s:2cosh(O0 = ( J - 1) ^Axi) + V ■ {e(F - 



+ (J- l)K2sinh(03), (5.36) 
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where the function ^ is defined by use of the mean value expansion sinh(0) = sinh(03) + 
cosh(^)(0 — 03). Unlike its counterparts in the analysis for the first two steps, this func- 
tion ^ is not piecewise smooth since 0^' of Eq.(1.4) belongs to >V^'^(i7) hence is only 
piecewise uniformly differentiable. The resulting mean value function ^ is therefore a 
piecewise uniformly continuous function. Because of the appearance of remaining sin- 
gular charges in the right hand side of Eq.(5.36), we know that is not in globally. 
Again we employ the decomposition = G/ + + 0*" to separate the singular compo- 

nents Gf,Gr and the regular component 0^. The first singular component Gj = Gjj 

3 

is induced by all Nj singular charges in 

-V-(e™FVG/) = {J (5-37) 

i 

Nr 

while the second singular component G^ = Grj is caused by all A^^ singular charges 

j 

in ^rnr 

Nr 

-V-(e„FVa) = (J-l)5^g,5(x,), (5.38) 

j 

and both singular components have estimates similar to Eqs.(4.1 1) and (4.12) 

11^ II ^ 11^- l\\Lo<'(n)NfKq^a^ 

^/ 

\\Gr\\L--in,) < j-^ in (5.40) 

||VG/i|L-(n,) < in (5.41) 

||VGr||L->(n,) < ^ in (5.42) 

By subtracting the singular components G/, Gr from Eq.(5.36) we obtain an equation for 
the regular component 

- V ■ (eFV0'^) + Jk^ cosh(O0" = V ■ (e(F - I)V03) + ( J - sinh(03) 

-V-((e-e^)FVG/) 
-V-((e-e^)FVa), (5.43) 

where the last two items V ■ ((e — em)FVG/) and V • ((e — em)FVGr) prescribe two 
interface conditions on the molecular surfaces Tj and F^: 

fcf = (e.-ejFVG/-n, (5.44) 
fcr = (e, - e„)FVa ■ n, (5.45) 
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similar to that defined in Eq.(4.28). For the regular component 0'', Theorem (4.4) states 
that it can be estimated with respect to the W^'^ norm as follows 



W\\w^,p{n) < C \^\\(l)''\\LP(n) + \\fGf\\w^-i/v.P{rf) + WfcrW VKi-i/p,p(r^.) 

+ ll/G/|lvyi-i/p,p(r^) + WfcrWw^-^/P'PiTr) 

+ IK J - l)/€^inh(03)||Lnn) + ||V • (e(F - I)V03)|| Lp{n) ) 

< C (^||/G/|lvKi-i/p.p(r/) + ll/Grlli4^i-i/p,p{r;) 

+ ll/G/||vyi-i/p,p(r^) + II/g,. ||vyi-i/p.p(rr) 

+ 1) sinh(03)||ip(Q) + ||F - I||i^i,p(n)||03|k2.p(Q)) 

< C (II J - l||vi/i.p(Q,)||F||iyi,p(n^) (5.46) 
+ II J - l||Loo(n)|| sinh(03)||ip(Q) 

+ l|F - I||H^i,p(n)||03||w2.p(n)) 
= C(|| J — l||^yi,p(Q^) + II J — l||Loo(f^^) + ||F — I||i4/i,p(Q)), (5.47) 

where in the last inequality we applied the estimate in Eq.(4.15) for the interface condi- 
tions fcj, and fcr- Finally, we estimate the change of the electrostatic forces due to the 
elastic deformation. By definition, the body force change is attributed to the variation 
of regular component(reaction field) 0'' and the variations of the the singular compo- 
nents(Coulomb potential field), and thus can be estimated as: 



ibsWLPin^f) 



< 



< C{\\J — l\\w^'P{ns) + ll-'^ — M\L°°{ns) + l|F — I||H/i.p(n))- (5.48) 
The surface force change in this step is defined to be 

fs-fs3 = -^e,(|FV0.p- |V03.r)n+^e^(|FV0™|2- |V03r„r)n 
— K^(cosh(0s) — cosh(03s)n 
= -^e,(FV0, -V03s)-(FV0, + V03s)n 

+ ^em(FV(^m - ^hm) ' (FV0m + V03m)n 

-/€2sinh(e')(0s-03s)n. (5.49) 

It follows that 

||fs — fs3|lH/i-i/p,p(rj) < c(^\\FV(f) — V</)3||H/i.p{n,) + ||FV</) — V(f)3\\wi'P(n'^) 

+ \\4>- (psWw^-pin,))- (5.50) 



To relate the estimate of FV0 — V03 to that of — 03 (the latter has already been 
estimated in Eq.(5.47)), we make use of the relation 

||FV0- V03||H/i.P(n.) = IIFV0-FV03 + FV03- V03iki.P(Q.) 

< ||F(V0 - V03)lki.p(n.) + ||(F - I)V03|ki.P(Q.) 

< l|F||iyi,P{n,)||(V0 - V03)||tyi,P(n,) + 
||(F - l)\\w^.p{n,)\\<p3\\w2'P{n,) 
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and a similar relation for ||FV0 — V03|| vi/i.p(o'j- By collecting these results together we 
can conclude from Eq.(5.50) that 

||fs — fs3||vyi-i/p.p(r^) < C^^ll-'^ — l||vi^i.p(n) + ||F — I||H/i>p(f7)j7 (5.51) 

which indicates the dependence and the boundedness of this electrostatic force compo- 
nent with respect to the elastic displacement field. 

5.5. Complete estimation of the electrostatic forces. The complete estimation of the 
electrostatic surface force is presented this lemma: 

Lemma 5.4. The electrostatic force can be made arbitrary small by reducing the vari- 
ations of ionic strength, the volume of the additional low dielectric space, the added 
singular charge and the magnitude of the elastic deformation. 

Proof. Following from its decomposition schemes (5.13,5.14), the estimation of total 
electrostatic body force and surface fore can be readily completed by combining their 
respective four components estimated in the four subsections above. The estimates for 
these two forces have an identical form 



||ff> — im\\LP(Q.^f) < — Kq\ + 1102 — 0l||ty2,p(Q^) + 1102 — <Pl\\w'^'P{n'^) 

+ Vmr + ||F — I||vyi>p(r2) + — l||vKi>p{n) j , (5.52) 

|fs — fso|lvKi-Vp,P(r^) < C{\^K — Kq\ + 1102 — 0l||ty2,p(Q^) + 1102 — 0l||w^2.P(f7^) 

+ Vmr + ||F — I||vFi.P(f7) + 11-^ ~ l||w^i>P(n) ) • (5.53) 



It is noticed in Eq.(5.33) that both ||02 — (j)i\\w'i,p(Q.^) and ||02 — 0i ||vF2.p(n;,) can be 
made arbitrarily small by adjusting the charges of added molecule Vtmr- Moreover 
F(0)(x) = 0, J(0)(x) = 1 follow from their definitions and both functions are infin- 
itely differentiable in the neighborhood of each function in 

Xp = {u G w^'^{nmf)\Mw^,P(n^,) < M}. 

Applying the Taylor inequality we have 

||F — I||iyi,P(0 < |||-DF|||||u||p^2,p(Q^) 

11-^— < II l-DJ| II vi/2,p(n3)! 

hence the last two items in estimates (5.52, 5.53) are also small for properly chosen 

Xp. □ □ 

6. Main Results: Existence of Solutions to the Coupled System 

We now establish the main existence result in the paper. It is noticed that for every 
element v G Xp one can derive a Piola transformation and solve for a unique potential 
solution of the Poisson-Boltzmann equation with this Piola transformation. The electro- 
static forces computed from this potential solution belongs to W^~^^^''p{T j) hence there 
is also a unique solution u to Eq.(3.1) corresponding to these electrostatic forces. This 
loop defines a map 5* which associates every v with a new displacement function u. Our 
existence result is based on the following version of the Schauder fixed-point theorem. 

Theorem 6.1. Let Xp be a closed convex set in a Banach space X and let S be a contin- 
uous mapping ofXp into itself such that the image of S{Xp) is relatively compact. Then 
S has a fixed-point in Xp. 
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Proof. See [38]. □ 

The Schauder Theorem depends on establishing continuity and compactness of the 
map S : Xp ^ Xp. We notice that Xp is convex and is weakly compact in W'^'^. 
Therefore the mapping S has at least one fixed-point in Xp if we can verify that S is 
continuous in some weak topology. 

Theorem 6.2. S : Xp Xp is weakly continuous in W'^''^{Qrnf)- 

Proof. This proof follows the similar arguments in [6]. Let v„ be a sequence in Xp 
and v„ ^ V in W^'^ as n — )■ oo. With these displacement fields, we can compute 
the electrostatic potential 0„ = 0(v„)( hence the electrostatic body force ffc„ = fb(vn) 
and surface force f^n = fs(vn)) and new displacement fields u„ = u(v„) defining the 
mapping S. We know from (5.52), (5.53) and (3.2) that 0„, f;,„, f^^ and u„ are bounded 
independently of n. Therefore there exists a subsequence v„j C v„, an electrostatic 
potential 0, and a displacement field u, such that 

0(v„J ^ as / — > CXD 

u(v„J ^ u as / — )■ oo 

We shall prove that u(v) = u by investigating the limit of the equations for 0(v„J and 
u(v„J, and of the expression for F(v„J, ffc(v„J, f,(v„J and J(v„J. Since $(v„J ^ 
$(v) in the same weak topology as v„j ^ v and W'^'^ is compactly embedded in C^, 
there is a subsequence of v„j. C v„j such that 

Vjifc — )■ V in C'^{il) as /c -> oo, 

•^(vnj ^i'v) in C^(fi) as A; oo; 

hence 

JiynJ J(v) in C^{Q) as A; oo, 

following the definition of J(v). The convergence of 

F(v„J ^ F(v) in C\n) as A; ^ oo, 

which involves the inversion of V$(v), is substantiated by continuous mapping from a 
n X n matrix to its inverse in C°(fi), i.e., 

in the neighborhood of each invertible matrix of C°(i7), and by the invertibility of 
V$(v„^) in W^'^{il). Now we can pass the equations satisfied by (f)n^ and v^^ to the 
limit and deduce that 

0(v„J ^ 0(v) = 0, 
u„j ^ u(v) = u. 

This proves the continuity of mapping S in the weak topology of W'^'^. □ □ 

Finally we verify that S{Xp) C Xp. By connecting the force estimates (5.52)-(5.53) 
and the estimate of displacement u in theorem (3. 1) we observe that 

||u||iy2,p < C'(^||ffe||LP(f7™./) + ||fs||vi/i-Vp,p(rj)^ 

< c(^\k- ko|||02 - (piWw^.piQs) + l|F - Illiyi.pcn,) + ll-'^ - j 
C 
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for appropriately small change in ionic strength and in the charges in the added molecules, 
where C, CI, C2 are the constants prescribed in inequality (3.5). Thus we verified that 
S{Xp) C Xp and $(u) is invertible. This gives the main result in the paper as the 
following theorem. 

Theorem 6.3. There exists a solution to the coupled nonlinear PDE system (3.1) and 
(1.4) for sufficiently small k — and sufficiently small rigid molecule Q^r with suffi- 
ciently small charges. 

Proof. This follows from Theorem 6.1 combined with Theorem 6.2. □ 

7. Variational Principle for Existence and/or Uniqueness 

In addition to the fixed point arguments, variational principles and quasivariational 
inequalities are also widely used for analyzing coupled systems of PDEs arising from 
multiphysics modeling. While quasivariational inequalities are exclusively used for sys- 
tems with boundary conditions given by inequalities, a single energy functional for the 
entire system is generally required for the application of either of these two approaches, 
and the stationary point of this energy functional with respect to each function shall pro- 
duce the corresponding differential equations and all boundary conditions. This energy 
functional is usually given by the total potential energy of the system, or by the sum of 
the potential energies of each equation if these energies are compatible. While it remains 
challenge to construct the total energy for our problem, we can give a coupled weak form 
of the entire system: 

{Ax,y) = (f6, v)i2(n,„^) + (f„v)z.2(r„^.) G P, (7.1) 

where x = (u, 0), y = (v, ip) are in the product space P of for the displace- 

ment field u and the W'^'P{fl) for the regular component of electrostatic potential 0*", i.e., 
P = VV^'Pinmf) X W^^P(^n), and the operator A is defined by 

{Ax,y) = (T(u),E(v)) + (eFV0,V0) + (Jfi:2sinh(0 + G'),V'), (7.2) 

where the stress tensor T and the strain tensor E were given in Eq. (1.3). 

Unlike the piezoelectric problems to which variational principles and quasivariational 
inequalities can be readily applied, we lack the coupling of the electrostatic potential and 
elastic displacement at the level of constitutive relations of the material [39, 41]. Instead, 
our electro-elastic coupling is through the electrostatic forces. We note that variational 
principles have been formulated for a class of fluid-solid interaction systems [42], which 
resemble our problem in that the coupling is through the boundary conditions of the 
elasticity equation instead of the constitutive relations. This will be examined for our 
problem in a future work. 

8. Concluding Remarks 

In this paper we have proposed and carefully analyzed a nonlinear elasticity model of 
deformation in macromolecules induced by electrostatic forces. This was accomplished 
by coupling the nonlinear Poisson-Boltzmann equation for the electrostatic potential field 
to the nonlinear elasticity equations for elastic deformation. The electrostatic of this cou- 
pled system is desribed by an implicit solvation model, and the Piola transformation de- 
fined by the solution of the elasticity equation is introduced into the Poisson-Boltzmann 
equation such that both equations can be analyzed together in a undeformed configura- 
tion. A key technical tool for coupling the two models is the use of an harmonic extension 
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of the elastic deformation field into the solvent region of the combined domain. Com- 
bining this technical tool with regularization techniques established in [2] and a stan- 
dard bootstrapping technique, we showed that the Piola-transformed Poisson-Boltzmann 
equation is also well-defined and the regular component of its solution has a piecewise 
VT^'^-regularity. This regularity matches that of the elastic deformation, giving access 
to a Schauder fixed-point theorem-based analysis framework for rigorously establishing 
the existence of solutions to this coupled nonlinear PDE system for small perturbation in 
the ionic strength and for small added charges. The existence of large deformation for 
large perturbations in ionic strength and/or charges can be obtained by combining our 
local result with general continuation techniques for nonlinear elastic deformation [36]. 
Our Shauder-type existence proof technique did not require that we establish a contrac- 
tion property for the fixed-point mapping S; this results in losing access to a uniqueness 
result for the coupled system, as well losing access to a fixed error reduction property for 
numerical methods based on the fixed-point mapping S. 

The coupling of elastic deformation to the electrostatic field is of great importance in 
modeling the conformational change in large macromolecules. To put this into perspec- 
tive, more comprehensive and realistic continuum models for macromolecular confor- 
mational changes can be developed based on the results in this article, for example, by 
coupling the (stochastic) hydrodynamical forces from the Stokes or Navier-Stokes equa- 
tion, or including van der Waals forces between closely positioned molecules. While 
mathematical models and robust numerical methods have been well studied for steady 
state fluid- structure interaction problems [6], the inclusion of van der Waals forces ap- 
pears to be more straightforward [34]. A major concern in applying these coupled mod- 
els, however, is the determination of the elasticity properties of macromolecules within 
the continuum framework, which requires new theoretical models and quantitative com- 
parisons between the continuum modeling and the classical molecular dynamical simu- 
lation and/or experiential measurements. In a future work we will study the development 
of numerical methods for this coupled system and apply this model to macromolecular 
systems where electrostatic forces play a dominant role. 
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